/*
 *  svk.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde on 1/04/2011.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "svk.h"

Svk::Svk()
:	HyperelasticInvariants()
{}

// ***

Svk::Svk(ifstream &in, const string& na, const string& t)
:	HyperelasticInvariants(na,t)
{
	string keyword;
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
		{
			if(!keyword.compare("E"))	in >> E;
			else if(!keyword.compare("NU"))	in >> nu;
	        else if(!keyword.compare("G"))	in >> mu;
	        else if(!keyword.compare("LAMBDA"))	in >> lambda;
	        else if(!keyword.compare("K"))	in >> K;
			else if(!keyword.compare("DENSITY")) in >> density;
		}

		if (E != 0.0 and nu != 0.0)
	    {
			mu = E / ( 2.0 * (1.0 + nu) );
	        lambda = nu * E / ((1.0 + nu)*(1.0 - 2.0 * nu));
	    }
	    else if (mu != 0.0 and nu !=0.0 )
	    {
	    	E = 2.0 * mu * (1.0 + nu);
	    	lambda = 2 * mu * nu / (1.0 - 2.0 *  nu);
	    }
	    else if (mu != 0.0 and lambda !=0.0 )
	    {
	    	E = mu * (3.0 * lambda + 2.0 * mu) / (lambda + mu);
	    	nu = lambda / 2.0 / (lambda + mu);
	    }
	    else if (E != 0.0 and K !=0.0 )
	    {
	    	mu = 3.0 * E * K / (9.0 * K - E);
	    	lambda = 3.0 * K * ( 3.0 * K - E) / (9.0 * K - E);
	    }
	    else if (mu != 0.0 and K !=0.0 )
	    {
	    	lambda = (3.0 * K - 2.0 * mu) / 3.0;
	    	E = 9.0 * K * mu / ( 3.0 * K + mu);
	    	nu = ( 3.0 * K - 2.0 * mu) / ( 6.0 * K + 2.0 * mu );
	    }
	    else if (E != 0.0 and mu != 0.0)
		{
			nu = (E - 2.0 * mu) / ( 2.0 * mu );
		    lambda = mu * (E - 2.0 * mu) / (3.0 * mu - E);
		}
		else if (E != 0.0 and lambda !=0.0 )
		{
		    mu = (E - 3.0 * lambda + sqrt(E*E + 9.0 * lambda * lambda + 2.0 * E * lambda)) / 4.0;
		    nu = 2.0 * lambda / (E + lambda + sqrt(E*E + 9.0 * lambda * lambda + 2.0 * E * lambda));
		}
		else if (nu != 0.0 and K !=0.0 )
		{
		    E = 3.0 * K * ( 1.0 - 2.0 * nu);
		    mu = 3.0 * K * ( 1.0 - 2.0 * nu ) / 2.0 / (1.0 + nu);
		    lambda = 3.0 * K * nu / ( 1.0 + nu);
		}
		else if (nu != 0.0 and lambda !=0.0 )
		{
		    mu = lambda * ( 1.0 - 2.0 * nu) / 2.0 / nu;
		    E = lambda * ( 1.0 + nu ) * ( 1.0 - 2.0 * nu ) / nu;
		}
		else if (K != 0.0 and lambda !=0.0 )
		{
		    mu = 2.0 / 3.0 * (K -lambda);
		    E = 9.0 * K * (K -lambda) / ( 3.0 * K - lambda);
		    nu = lambda / (3.0 * K - lambda);
		}

	    if (K == 0.0)
	    {
	        K = lambda + 2.0/3.0 * mu;
	    }

	    cout << "Saint-Venant Kirchhoff Material constants have been set to:" << endl;
	    cout << "\t" << "E= " << E << endl;
	    cout << "\t" << "nu= " << nu << endl;
	    cout << "\t" << "mu= " << mu << endl;
	    cout << "\t" << "lambda= " << lambda << endl;
	    cout << "\t" << "K= " << K << endl;
	
	Ct(0,0) = lambda + 2.0 * mu;
	Ct(0,1) = lambda;
	Ct(0,2) = Ct(0,1);
	
	Ct(1,0) = Ct(0,1);
	Ct(1,1) = Ct(0,0);
	Ct(1,2) = Ct(0,1);
	
	Ct(2,0) = Ct(0,1);
	Ct(2,1) = Ct(1,2);
	Ct(2,2) = Ct(0,0);
	
	Ct(3,3) = mu;
	Ct(4,4) = Ct(3,3);
	Ct(5,5) = Ct(3,3);
	
	hypertype = "coupled";
}

// ***

void Svk::strainEnergy()
{
      if (J > 0.0)
      {
    	  W = lambda / 8.0 * (I - 3.0) * (I - 3.0) + mu / 4.0 * (I * I - 2.0 * I  - 2.0 * II + 3.0);
      }
      else
      {
    	  W = -1.0;
      }
}

// ***

void Svk::DstrainEnergy()
{
      if(J > 0.0)
      {
    	  W1 = lambda / 4.0 * (I - 3.0) + mu / 2.0 * (I - 1.0);
    	  W2 = - mu / 2.0;
    	  W3 = 0.0;
      }
      else
      {
    	  W1 = -1.0;
    	  W2 = -1.0;
    	  W3 = -1.0;
      }
}

//***

void Svk::DDstrainEnergy()
{
      if(J > 0.0)
      {
    	  W11 = lambda / 4.0 + mu / 2.0;
    	  W12 = 0.0;
    	  W13 = 0.0;
    	  W22 = 0.0;
    	  W23 = 0.0;
    	  W33 = 0.0;
      }
      else
      {
    	  W11 = -1.0;
    	  W12 = -1.0;
    	  W13 = -1.0;
    	  W22 = -1.0;
    	  W23 = -1.0;
    	  W33 = -1.0;
      }
}

//***

void Svk::DDDstrainEnergy()
{
    if(J > 0.0)
    {
        W11 = lambda / 4.0 + mu / 2.0;
        W12 = 0.0;
        W13 = 0.0;
        W22 = 0.0;
        W23 = 0.0;
        W33 = 0.0;
    }
    else
    {
        W11 = -1.0;
        W12 = -1.0;
        W13 = -1.0;
        W22 = -1.0;
        W23 = -1.0;
        W33 = -1.0;
    }
}

//***

//double Svk::getE()
//{	return E;}
//
//// ***
//
//double Svk::getNu()
//{	return nu;}
//
//// ***
//
//double Svk::getLambda()
//{	return lambda;}
//
//// ***
//
//double Svk::getMu()
//{	return mu;}
//
//// ***
//
//double Svk::getDensity()
//{	return density;}
